Novel connections and physical implications of thermal metamaterials with imperfect interfaces

Thermal metamaterials are of great importance in advanced energy control and management. Previous studies mainly focused on interfaces with perfect bonding conditions. In principle, imperfectness always exists across interface and the effect is intriguingly important with small-length scales. This work reports the imperfect interface effect in thermal metamaterials thoroughly. Low conductivity- and high conductivity-type interfaces are considered. We show that an object can always be made thermally invisible, with the effect of imperfect interface, as that of a homogeneous background material. This unprecedented condition is derived in an exact and analytic form, systematically structured, with much versatile and physical implications. Conditions for thermal shielding and enhancements are analytically found and numerically exemplified, highlighting the specific role of material and geometric parameters. We find that both types of interfaces are complementing with each other which, all together, will constitute a full spectrum to achieve the thermal invisibility. The analytic finding offers a general perception that adds to the understanding of heat transport mechanism across interfaces in thermal metamaterials, in ways that drastically distinct from that of ideal interfaces. This finding opens up new possibilities for the control and management of thermal metamaterials with imperfect bonding interfaces.

Management and control of heat is familiar to all of us and is ubiquitous in our daily life. The phenomena of heat transfer across solid interfaces between different materials are of central importance in modern technologies and applications, particularly with scaling down of the devices 1 . On the other hand, thermal devices that are capable of confine energy, preferably with an enhanced intensity, are becoming imperative especially in the field of energy harvesting. To realize these functionalities in various applications, it would require an in-depth understanding on the transfer mechanism across interfaces 2 , particularly counting into the effect of bonding imperfectness.
In the last decade, significant progress was made in the compositions of materials, referred to as thermal metamaterials, to control the flow of heat along a certain path. Among them, by surrounding an artificially designed material, it is possible to make an object undetectable from thermal measurements 3,4 . The coordinate transformation method 5 , originally proposed in optics for invisibility cloaks, offers a systematic approach capable of precisely manipulating wave-based as well as diffusion-based phenomena 6 in a desired manner. Han et al. 7 showed that a circular anisotropic layer can serve as a tunable device to control thermal energy within a certain area. Other relevant works on thermal devices have been studied theoretically and numerically [8][9][10][11][12] . Substantial progress has also been made in experimental verifications of thermal devices [13][14][15][16][17][18][19][20][21][22][23] . For a general exposition of the subject, one can refer to the review articles [24][25][26][27][28] and the monograph 29 . All these aforementioned studies mainly consider ideal interfaces. In principle, imperfectness of bonding always exists across the interface between two dissimilar solids in contact. Particularly, in small-length scales, an ideal interface framework could not be adequate and that interfacial thermal contact will become non-negligible 2,30 . For a given volume fraction, the smaller the partcle size, the larger will be the total contact surface area between different materials. Yet very little work focused on the effects of imperfect interface effect, nor on the general transport mechanisms of thermal energy across interfaces of thermal metamaterials. Zheng and Li 31 recently demonstrated that, with the interfacial thermal resistance, the temperature field outside the cloak could be greatly distorted, in contrast to that of ideal interfaces.
In this study, we consider steady-state heat conduction, under a macroscopic scale, governed by Fourier's law. We examine the effect of imperfect interface in the design of a homogeneous thermal metamaterial, mostly

Results
At a macroscale, heat conduction can be viewed as a diffusion process in which thermal energy flows from hot region to cold region. The heat equation, described by Fourier's law and the conservation of energy, characterizes the temperature field over space and time. Under a steady state condition, in the absence of heat source, the temperature field in a medium is governed by the Fourier law ∇ · (k∇T) = 0 , where T is the temperature field and k is the thermal conductivity tensor. For natural materials the values of thermal conductivity vary a great deal, ranging from poorly conducting materials with k ≈ 0.02 W m −1 K −1 to good conductors with k ≈ 400 W m −1 K −1 .
We consider a two-dimensional model depicted as in Fig. 1. Region I is isotropic with conductivity k = k 0 , which is the region for heat flux control. Region II is the thermal metamaterial made of an anisotropic layer with constant radial and circumferential conductivities given as Region III is the surrounding background material with the same isotropic conductivity k = k 0 as in Region I. Along the x-direction, a uniform thermal intensity E is prescribed. The interfaces between Region I and Region II, and between Region II and Region III, could be imperfectly bonded, with either LC-or HC-type. We will show that a tailored set of constant materials and geometric parameters in Region II can act as a functional thermal metamaterial in such a way that the outer field in Region III will remain to have a constant temperature intensity E, as if the whole medium is spatially homogeneous. In addition, the thermal intensity in Region I can be tuned to be greatly enhanced or shielded.
Low conductivity-type interface. The LC-type interface is that the interface exhibits thermal contact resistance, or Kapitza resistance 30 , between the two surfaces in contact. In this case, the temperature will have a drop across the interface, while the normal component of heat flux q · n = −(k∇T) · n will be continuous. For theoretical and numerical studies in the context of composites, we mention the works of Refs. [39][40][41][42][43] . Specifically, the interface conditions for the LC type of interface at r = a are www.nature.com/scientificreports/ The interface condition at r = b will be given in Supplementary Material S3 for brevity. Here ∂/∂n is the normal derivative at the interfaces and T 1 , T 2 and T 3 represent the temperature fields in Regions I, II and III, respectively. The interface parameter β is defined as the ratio between the temperature drop and the heat flow across the interface, with a dimension of W m −2 K −1 . When β → ∞ , the LC-type interface will reduce to that of the perfect bonding interface. We shall use β a to denote the LC-type interface parameter at interface r = a, and similarly β b for r = b. It should be emphasized that the mathematical framework of LC-type interface can be formulated by considering a thin interphase layer, of a constant thickness t with relatively low conductivity k c , compared to those of the two neighboring regions 39 . By taking an asymptotic analysis with equilibrium considerations, one can show that the effect of the interphase is equivalent to an interface between the two phases with a LC-type of interface condition, where β = lim t→0,k c →0 k c /t . We will show that, in the presence of thermal contact resistances, it is possible to achieve the thermal invisibility in the outer Region III. The exact connection for the LC-type can be exactly determined as where c = (a/b) 2 is the area fraction of Region I, g = k 0 /k G , β a = k 0 (aβ a ) and β b = k 0 (bβ b ). The quantity k G = √ k r k θ is the geometric average of k r and k θ , and = √ k θ /k r is a parameter that signifies the anisotropy of the conductivity tensor in Region II. Apparently, λ = 1 stands for an isotropic material in Region II. We note that all the parameters in (3) are non-dimensional. Here the top tilde of β is used for the non-dimensional interface parameter with LC-type interfaces, while a top hat will be used for those of HC-type interfaces to be described later. The constant g is a parameter indicating the relative magnitude of k G and k 0 . For example, g < 1 means k 0 < k G , whereas g > 1 indicates k 0 > k G . Note that g = 1 is exactly the thermal neutrality condition for the ideal interfaces. In our formulation, we consider that both the interfaces at r = a and r = b are imperfectly bonded with different interface parameters. When only one interface is imperfect, we can simply let β a = 0 or β b = 0 in (3) and directly deduce a simplified version of (3). Particularly, when the thermal invisibility condition (3) is fulfilled, the temperature field in the Region III will remain as T 3 = −Er cos θ . This suggests that Regions I and II, togther with the effect of the imperfect interfaces, are effectively equivalent to that of a homogenous circular cylinder with conductivity k 0 . We note that, in addition to the interface parameters β a and β b , the thermal invisibility condition (3) also depend on three parameters, g = k 0 /k G , λ and c = (a/b) 2 . Note that in the case of perfect bonding interfaces, we have β a → ∞,β b → ∞ , that is β a → 0 , β b → 0 . Equation (3) will exactly reduce to the effective conductivity of composte cylinder assemblage 44 , an isotropic inclusion of conductivity k 0 coated with a cylindrically orthotropic layer with conductivity tensor (1). It can be verified analytically, by direct expansion, that the thermal invisibility condition (3) will recover exactly to the known neutrality condition for the perfect bonding situation 7 , k rr k θθ = k 2 0 . Note that this latter connection under the perfect bonding condition is independent of the absolute sizes of the regions I and II, nor of the area fraction. Physically, the LC-type interface, compared to that of perfect interface, exhibits a thermal contact resistance at interface. Thereby, in the presence of a positive interface parameter β, to fulfill the thermal invisibility condition, one would necessarily have k G > k 0 , i.e., g < 1.
High conductivity-type interface. The behavior of HC-type interface is physically opposite to that of LC-type interface. The mathematical framework of HC-type interface can be derived by considering a thin interphase layer of a constant thickness t with relatively high conductivity k c . Related studies on the effect of HC- Figure 1. A schematic illustration of a circular thermal device with imperfect interfaces. Region I has an isotropic conductivity k = k 0 and is the targeted region for heat flux control. Region II is the designed region with constant radial and circumferential conductivities. Region III is the surrounding background material with the same isotropic conductivity as in Region I. The interfaces at r = a and r = b are imperfectly bonded with different interface parameters. www.nature.com/scientificreports/ type interfaces in composite systems can be found, for example, in the papers 33,[45][46][47][48][49] . In the context of elasticity, HC-type interface is relevant to the interface conditions with surface stress effects 50 . Specifically, the HC-type interface conditiond at r = a are where s is the surface Laplacian operator. The interface condition at r = b is omitted here and is expressed in (S4) for conciseness. The parameter of the HC-type interface can be derived as α = lim t→0,k c →∞ k c t , where α has a dimension of W K −1 . We note that α a = α b = 0 will correspond to ideal interface conditions. For the configuration of Fig. 1, we find that, in the presence of HC-type interfaces, when the material and geometric parameters fulfill the condition the temperature field in Region III will remain to have T 3 = −Er cos θ , as if the whole medium is homogeneous.
In (5), we have defined the dimensionless parameters α a = α a (ak 0 ) , α b = α b (bk 0 ) . Note again, as in (3), in addition to the interface parameters α a and α b , the invisibility condition (5) also depend on the three parameters, g, λ and the area fraction c. Note also that for perfect bonding interfaces we have α a → 0 , α b → 0 , namely α a → 0 and α b → 0. It can be shown that Eq. (5) will recover the known condition under the perfect bonding conditions 7 k rr k θθ = k 2 0 . Physically, the HC-type interface represents a more conductive behavior than that of the perfect bonding interface. As such, when we have k 0 > k G , to fulfill the invisibility condition it is necessary to have HC-type interfaces to balance the mismatch between k 0 and k G . To address the point, in Fig. 2 we plot the relation between k 0 /k G and the dimensionless interface parameters, β a or α a . In the numerical illustrations, we have assumed β a = β b = β and α a = α b = α, as both interfaces separate Region II and a region with isotropic conductivity k 0 . For numerical illustrations, we specify c = (a/b) 2 = 1/4. Note that the horizontal axis is expressed in logarithmic scale. The left panel of Fig. 2 corresponds to the range of k 0 /k G < 10 0 , where a LC-type interface is needed to compensate the conductivity of Region II so that the outer field in Region III will remain unaffected.
The associated dimensionless interface parameter k 0 /aβ is indicated on the vertical axis on the left panel. On the other hand, when k 0 /k G > 10 0 , a HC-type interface will be invoked to render the outer field invisible, as that of a homogeneous medium. The corresponding non-dimensional interface parameter α/ak 0 is labelled on the vertical axis on the right side. Different curves for various values of λ are plotted. The curves corresponding to λ = 10 and λ = 10 2 are nearly overlapping, and thus for λ > 10 2 the corresponding curve can be viewed as an asymptotic limit. We note that all curves are intersecting at a single point k G /k 0 = 10 0 , which is the situation for www.nature.com/scientificreports/ perfect bonding interfaces corresponding to α → 0 or β → ∞. Interestingly, the curves between the two sides along the line k 0 /k G = 10 0 appear to possess a reflection symmetry. We will confirm this later in "Methods" section. In Fig. 3, we plot the value of β b = k 0 bβ b versus the anisotropy parameter λ for the case of g = 2/3, c = 1/4 and b = 1 μm. It is seen that when λ → 0, we have β b → 0 . On the other hand, when λ → ∞, we have β b → 1 3 . Also we observe when λ < 10 −2 , or λ > 10 −2 , the curves are virtually unchanged. Note Fig. 3 also applies to the case of g = 3/2, but now the vertical axis needs to be replaced by the scale of α b . In short, the exact thermal invisibility conditions are expressed in analytical forms in Eq. (3) for the LC-type and in Eq. (5) for the HC-type interfaces. The explicit expressions depend on the material properties as well as the geometric parameters. Note that all these parameters are non-dimensional, including the interface parameters β and α , the anisotropic ratio = k θ k r , the ratio g = k 0 /k G and the area fraction c = (a/b) 2 . A detailed study on the effects of extreme values of , g = k 0 /k G and c = (a/b) 2 is given in Supplementary Materials. This will provide theoretical examinations, in complement with the numerical illustrations given in Figs. 2 and 3, on the trend under the extreme values of g = k 0 /k G , and c to achieve the thermal neutrality conditions (3) or (5). Both theoretical and numerical results agree perfectly under extreme values of g for Fig. 2 and for Fig. 3.

Theoretical analysis
We now present the theoretical analysis of the conduction of the geometric model shown in Fig. 1 ∂θ 2 = 0. Boundary condition T(x = r cos θ = ±x 0 ) = ∓T 0 are prescribed on the left and right surfaces of the configuration (Fig. 1), with T 0 = −Ex 0 and E is the prescribed constant intensity. Taking into account the boundary condition and the geometric symmetry, the temperature potentials in the three regions are given in (S2). If the medium is homogeneous throughout, the temperature gradient will be spatially uniform. But in the presence of the Region II together with the imperfect interfaces along r = a and r = b, heat flows will be generally distorted. We find that under the condition of (3) or (5), the outer field in Region III will remain to possess a uniform intensity E, as if the Region II and the imperfect interfaces are entirely invisible. A direct outline of the solution procedure will be given in Supplementary Material. In "Methods" section, instead, a simple exposition of the thermal invisibility condition based on the neutral inclusion concept is presented through a series of homogenization process (see Fig. 6). Now under the invisibility condition (3), we find that the temperature field in Regions I and II can be exactly expressed as where T 3 = −Er cos θ . Here the subscript β on the left-hand side is used to signify the fields for LC-type interfaces, in distinction with that of HC-type interfaces to be presented later. We are particularly interested with the ratio of T 1 /T 3 , in which a value of T 1 /T 3 > 1 will represent a concentrating effect, while a value of T 1 /T 3 < 1 will indicate a cloaking or shielding effect. We first note that when c → 0, that is a ≪ b , we have T 1 /T 3 → 0. On the other hand, when c → 1, that is a ≈ b, we have T 1 /T 3 → 1. Since 0 < c ≤ 1, it is evident that −∞ < ln c ≤ 0 and 1 ≤ c −1 / 2 = b a < ∞ . Also, from elementar y algebra, we know that cosh 2 ln c ≥ 1 and (6) www.nature.com/scientificreports/ −1 < tanh 2 ln c < 0 . Thus, in view of (6), we conclude that the constant intensity in Region I is always bounded by the ratio b/a, where the upper limit is attained by the perfect bonding conditions, that is g = 1,β a =β b = 0 , and λ → 0. Note that for perfect bonding interfaces, the internal fields will be exactly reduced to the known solutions To illustrate the variation of intensity in Region I, we see that the constant intensity in (6) will depend on c, λ, and 1 +β a g . Note these three parameters need to follow the thermal invisibility condition (3). For convenience, let us define � = 1 +β a g and rewrite (3) as in Eq. (S5). In Fig. 4, we plot the value of T 1 /T 3 versus λ for different values of Λ, where the area fraction c is specified as c = 1/4 . Since β a ≥ 0 and 0 < g < 1, the value of Λ is always greater than 1 and that Λ = 1 corresponds to the perfect bonding condition. If one defines cr as the value of λ such that T 1 /T 3 = 1, it can be seen from Fig. 4 that, for a fixed value of c, when the value of Λ is increased, a smaller value of λ is needed to attain the value T 1 /T 3 = 1. Next, we consider the HC-type interfaces. Now under the condition (5), we find that the temperature field in the Regions I and II can be exactly expressed as We note that the two terms in the denominator, as in (6) and (7), are always greater than unity, thus it is seen that T 1 /T 3 ≤ b/a. The upper limit of is achieved for λ → 0 under the perfect bonding condition, that is g = 1 and α a =α b = 0 . Again, the expressions (9) and (10) can be shown to reduce to the known solution for perfect bonding interfaces. For simplicity, we can define Ŵ = g 1 +α a and rewrite the thermal invisibility condition (5) as in (S6). The value of Ŵ is always greater than 1, since α a ≥ 0 and g > 1, and that Γ = 1 corresponds to the perfect bonding situation. For a fixed value of c, one can draw diagrams for the value of T 1 T 3 α versus λ for different values of Ŵ. Note that the illustration for HC-type interfaces will be exactly the same with that of the LC-type interface simply by replacing Λ by Γ, which will be explained later in "Methods" section.
In Fig. 5a,b, we plot the temperature and heat flux profiles based on the analytical solutions, (6), (7), (9), (10), and numerical simulations based on finite element calculations (COMSOL). The left panel in Fig. 5 is the analytical solutions based on (6) and (7) for the LC-type interfaces, and (9) and (10)for the HC-type interfaces. The right panel is numerical simulations based on finite element calculations (COMSOL). We mention that the finite element software COMSOL considers that the interfaces are perfect. To simulate the effect of imperfect bonding, we model the imperfect interface as a thin interphase layer of a constant thickness t with isotropic conductivity k C , prescribed by either k 0 t aβ a (LC-type) or α a ak 0 /t (HC-type). Here we consider a thin thickness t = b/10 3 = 10 −3 μm. In the numerical illustrations, we consider c = 1/4, b = 1 μm and the value of λ = 1. For the LC-type case, we consider g = 2/3, while for the HC-type interface we consider g = 3/2. To perceive the temperature cosh ln r a 1 + g 1 +α a tanh ln r a cosh 2 ln c 1 − g 1 +α a tanh 2 ln c .  (6) and (7) for LC-type interfaces, and (9) and (10) for HC-type interfaces. The right panel is the numerical simulations based on finite element calculations (COMSOL). In numerical calculations, we consider c = 1/4 and b = 1 μm. For g = 2/3, LC-type interfaces are invoked. Temperature contours for λ = 10 0 is illustrated in (a). For g = 3/2, HC-type interfaces are necessary, and contour plots for λ = 10 0 are given in (b). A temperature or heat flux profile along the x-axis is shown below the contour diagram to highlight the jump of relevant quantity across the interface. www.nature.com/scientificreports/ drop or heat flux drop across the interface, a profile along the x-axis is highlighted just below the contour diagram. It is seen that both solutions, analytical solutions and finite-element simulations, agree well. For completeness, additional demonstrations for the temperature contours are presented for λ = 10 −1 and λ = 10 in Supplementary Materials, which demonstrate the effect of material anisotropy in Region II. It is seen that the case of λ = 10 −1 will lead to a concentrating effect, while λ = 10 will result into a thermal shielding effect. In Supplementary Material Table S1, we also list a comparison of temperature field calculated from finite element numerical simulations, based on COMSOL with three different interphase thickness, t = b/1200, t = b/1000 and t = b/800, and the analytic results for the LC-type interfaces, (6), (7), or for the HC-type interfaces, (9), (10). The percentage error is estimated by (COMSOL prediction − analytic result)/analytic result. We see that the error percentage for t = b/1200 and t = b/1000 are nearly the same within 0.3% at most.

Methods
Thermal invisibility condition. The thermal invisibility condition (3) for the LC-type interfaces and (5) for the HC-type interfaces can be constructed via the concept of neutral inclusions 32 through a series of homogenizations via composite cylinder assemblages (CCA) 51 . A schematic illustration is presented in Fig. 6, in which the derivation can be separated into three steps. First, we consider a circular cylinder of radius a (Region I) with a LC-or HC-type of boundary surface. The theoretical formulation is to identify an effective conductivity for Region I, without the surface effect, so that the heat flow pattern outside the circular cylinder will be identical for both configurations. The effective conductivity can be exactly derived as k 0 = k −1 0 + (aβ a ) −1 −1 = k 0 1 +β a for a LC-type interface, and k 0 = k 0 + α a a = k 0 1 +α a for a HC-type interface. Next we consider a circular cylinder with conductivity of either k 0 or k 0 , coated with a circular layer of cylindrically orthotropic material with conductivity tensor (1). This is an exactly solvable composite cylinder assemblage that is amenable to an exact determination for the effective conductivity 51 . www.nature.com/scientificreports/ The effective conductivity of composite cylinder assemblages consisting of an isotropic core, with conductivity k 0 or k 0 , coated with a cylindrically orthotropic material with conductivity tensor (1) can be expressed as k CCA or k CCA 44 , where g =k 0 k G and ĝ =k 0 k G . The last step is to consider a circular cylinder with isotropic conductivity k CCA ork CCA incorporating with the interface effect with surface parameter of β b or α b . The analytical result is similar to that in the first step. But now, to be fully compatible to the surrounding Region III, the effective conductivity of the last step is necessarily identical to the surrounding medium, Region III, with conductivity k 0 .
A duality relationship. For a two-dimensional medium, a duality relation is known to exist in the context of conduction phenomena, which links the effective conductivity of a composite to that of a composite with the same phase geometry but replacing its phase moduli with the inverse conductivities 32,52-54 . The existence of the connection between the dual systems relies on the fact that a two-dimensional divergence free field, when rotated pointwise by a 90 degree, produces a curl free field, and vice versa. Lipton 45 proved that the duality relation is also applicable to a composite system with imperfect interfaces. Here we demonstrate that our invisibility conditions (3) and (5) can indeed be linked with each other via the duality relation. For simplicity of expositions, we will make use of Proposition 2 given in the work 55 , recorded as k * β (k 1 , k 2 , . . . , k n , β) = 1 k * α k −1 1 , k −1 2 , . . . , k −1 n , α −1 . Here the superscript k * β represents the effective conductivity of the composite system with LC-type interfaces, while k * α is the effective conductivity of the dual system with HC-type interfaces. The arguments inside the parenthesis of k * α or k * β are simply the phase conductivities and the interface parameters of the composite system. Now letting k 0 → k −1 0 , k −1 CCA →k CCA and β −1 → α in (12) 1 , we can see that (12) 2 is readily deduced. Likewise, one can deduce (12) 1 from (12) 2 . We thus conclude that the system with LC-and HC-type interfaces are complementary with each other, that provides a link to achieve thermal invisibility for the case of k G > k 0 or k G < k 0 .

Conclusions
Utilizing the theory of neutral inclusions in mechanics of composites, we provide a thorough analysis on the effect of imperfect interfaces in functional thermal metamaterials, made of homogeneous anisotropic conductivity. When the interfaces are perfectly bonded, the condition of k 2 0 = k r k θ will ensure that the temperature field in Region III will be the same as that of a homogeneous background material. While for interfaces, the neutrality conditions can be constructed in the well-structured forms (3) and (5). Specifically, they depend on the absolute size, the material anisotropy and the area fraction c. This finding now works at new levels of applicability, irrespective of the bonding perfectness g . In general, the imperfect interface will broadly affect the heat transfer across interface and will exhibit size effect. We demonstrate quantitatively the interplay between the material and geometric parameters of the configuration, and also explore the extent to which a thermal confinement or thermal shielding can be actually achieved within a targeted region. Numerical simulations based on finite element calculations are also carried out to validate our analytical results. We prove that the thermal invisibility conditions with the LC-and HC-type interfaces are reciprocal to each other which, all together, will span a full spectrum for the feasibility of thermal invisibility. We believe that the study provides further insights to the understanding of heat transfer mechanisms across interfaces in the performances of thermal metamaterials, and will motivate future experiments for practical realization. Since the interface in general may not be always perfect, this finding can serve as a general guideline in the design of thermal metamaterials that facilitates broad applications. (11) k CCA =k G g + 1 + c g − 1 g + 1 − c g − 1 ;k CCA =k G ĝ + 1 + c ĝ − 1 ĝ + 1 − c ĝ − 1 ,

Figure 6.
Homogenizations of the coated inclusion with imperfect interfaces. The first step is to consider Region I with the effect of imperfect interface at r = a. This will give a circular region with an effective isotropic conductivity k 0 or k 0 . Secondly, the composite cylinder model of an isotropic core with conductivity k 0 or k 0 coated with a cylindrically orthotropic material is determined, expressed as k CCA or k CCA . The final step is to homogenize a uniform conductivity k CCA or k CCA with an imperfect effect at r = b. This will lead to the thermal invisibility condition (3) and (5), respectively.